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Abstract 

A turbulence module is developed for the 2D ver- 
sion of the NPARC code which is currently restricted 
to planar or axisymmetric flows without swirling. 
Four turbulence models have been built into the 
module: Baldwin-Lomax, Chien, Shih-Lumley and 
CMOTT models. The first is a mixing-length eddy- 
viscosity model which is mainly used for initializa- 
tion of computational fields and the last three are the 
low Reynolds number two- equation models. Unlike 
Chien’s model, both the Shih-Lumley and CMOTT 
models do not involve the dimensionless wall dis- 
tance 3 /+, an advantage for separated flow calcu- 
lations. Contrary to the NPARC and most other 
compressible codes, the non-delta form of transport 
equations is used which leads to a simpler lineariza- 
tion and is more effective than using the delta form 
in ensuring the positiveness of the turbulent kinetic 
energy and its dissipation rate. To reduce numeri- 
cal diffusion while maintaining necessary stability, a 
second-order accurate and bounded scheme is used 
for the convective terms of the turbulent Jr ansport 
equations. This scheme is implemented in a deferred 
correction manner so that the main coefficients of 
the resulting difference equations are always posi- 
tive, thus making the numerical solution process un- 
conditionally stable. The system of equations are 
solved via a decoupled method and by the alternat- 
ing direction TDMA of Thomas. The module can 
be easily linked to the NPARC code for turbulent 
flow calculations. 

1. Introduction 

It has been long recognized that there is a 
gap between turbulence model developers and CFD 
users. The former mainly use simple flows to ver- 
ify new modeling concepts and evaluate the result- 
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ing models, while the latter are usually reluctant to 
implement and test new, more advanced turbulence 
models until they see such models showing good per- 
formance for a wide range of complex flow situations. 

In order to bridge this gap, we have developed 
turbulence modules for CFD codes used in industry. 
Under the widely used Boussinesq’s isotropic eddy- 
viscosity concept, the Reynolds-averaged equations 
governing turbulent flows are of the same form as 
those governing laminar flows, except that the lam- 
inar viscosity ft is replaced by the effective viscosity 

Me//=M + Mt (1) 

Therefore, a mean flow solver can be used to cal- 
culate turbulent flows once the turbulent eddy- 
viscosity fi t is available. The module is written in a 
self-contained manner so that the user can use any 
turbulence model in the module without concern as 
to how it is implemented and solved. The input to 
the module are the mean flow variables, boundary 
and geometric information which are to be provided 
by a mean flow solver. The output of the module 
are the turbulent eddy- viscosity and relevant tur- 

bulent source terms which are needed for the mean 
flow calculation. Here the turbulent source terms 
exist only for more sophisticated turbulence models 
beyond the level of isotropic eddy- viscosity models. 
For most CFD codes, especially those for compress- 
ible flow calculations, the laminar viscosity fi is a 
variable not a constant. In this case, few changes are 
required for a mean flow solver to use the turbulence 
module. The interaction between the mean flow 
solver and the turbulence module will give the final 
turbulent flow solution, as shown in Fig.l. With the 
aid of the modules, turbulence model developers can 
also take the advantage of the well-established and 
sophisticated CFD codes to test turbulence models 
for a variety of complex flows which are intractable 
with simple research codes. 

The NPARC code has been extensively used in 
the U.S. aerospace and aeropropulsion communities. 
The code can handle very complicated geometries. 
However until recently, only algebraic models such as 
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the Thomas and Baldwin-Lomax models have been 
available in the NPARC code for turbulent flow sim- 
ulations. The Chien’s low Reynolds number K-e 
model 1 with modifications for compressibility has 
been available in the 2D version but was not suc- 
cessfully installed in the 3D version. Another K-e 
model (the NPARC 1.0 K-e model) was installed 
in both the 2D and 3D versons but has not pro- 
vided desirable accuracy and stability. Recently, the 
two-equation turbulence model in the NAPRC code 
(both the 2D and 3D versions) has been modified so 
that the model is based on the low Reynolds number 
K-e model of Chien and no longer on the NPARC 
1.0 K-e model. Stability enhancements and a new 
inflow boundary condition for the turbulent quanti- 
ties were also added to the K-e model. Comparisons 
of the NPARC solutions obtained using the previous 
and newly installed models with experimental data 
indicated that the Chien’s K-e model installed im- 
proves the capability of the NPARC code for propul- 
sion turbulent flow calculations 2 . 

In the module developed for the NPARC code, 
the three low Reynolds number K-e turbulence mod- 
els have been implemented: Chien 1 , Shih-Lumley 3 
and CMOTT models. The Baldwin-Lomax mixing- 
length model is also included in the module, but it is 
mainly for initializing computational fields. Chien’s 
model is one of the well-known low Reynolds num- 
ber K-e models. However, it has some undesirable 
deficiencies. First, a near wall pseudo-dissipation 
rate is introduced to remove the singularity in the 
dissipation rate equation at the wall. The defini- 
tion of the near wall pseudo-dissipation rate is quite 
arbitrary. Second, the model constants are differ- 
ent from those of the standard K-e model 4 , making 
the near wall model less capable of handling flows 
containing both high Reynolds number turbulence 
and near wall turbulence, which is often the case 
for practical flows. Patel et al. 5 require as the first 
criterion, the ability of the near wall models to be 
able to predict turbulent free shear flows. Third, the 
dimensionless wall distance t/ + is used in the damp- 
ing function / M for the eddy- viscosity. Because the 
y + involves the friction velocity u T which is equal to 
zero at separation or reattachment points, any model 
using y + may have difficulties for separated flows. 
The Shih-Lumley and CMOTT models are free of 
these deficiencies. The two models differ from one 
another in the C ^ formulation, one using the stan- 
dard constant and the other a new formulation. The 
new C M has the following desirable features: (a) It 
is derived from a rigorous realizability analysis 6 that 
requires the non-negativity of the turbulent normal 
stresses and Schwarz’ inequality between any fluctu- 


ating quantities. As a result, unlike most of the ex- 
isting models, it satisfies the realizability conditions, 
(b) It accounts for the effect of the mean deformation 
rate by which the eddy- viscosity will be significantly 
reduced to an adequate level to mimic complex flow 
structures, (c) It is easier to use, as compared with 
other formulations. Simplicity is of great value for 
practical engineering applications. Successful appli- 
cations of this new can be found in Shih et al. 6 ’ 7 . 

Turbulence model equations require special 
treatment to ensure numerical realizability such as 
the positiveness of the turbulent kinetic energy and 
its dissipation rate. They we also often of source- 
dominate nature, which makes the linearization of 
source terms crucial for computational stability. 
Due to these considerations, we use the non-delta 
form of the transport equations in the module, which 
is contrary to the NPARC code as well as most other 
compressible codes. The non-delta form leads to 
simpler linearization and is more effective to ensure 
the positiveness of the turbulent kinetic energy and 
its dissipation rate than the delta form. For sim- 
plicity and also due to the fact that the coupling 
between the turbulence model equations is, in gen- 
eral, not very critical for the overall solution process, 
these equations are solved in a decoupled manner. 
The discretization of the convection terms is another 
important issue. The key point is to use a scheme 
which reduces numerical or false diffusion to such a 
level that it will not obscure the real viscous process, 
while maintaining necessary stability. Here, a hybrid 
linear/parabolic approximation (HLPA) scheme of 
second-order accuracy® is used. It is shown 9 that 
the HLPA scheme is capable of yielding low diffu- 
sive and always bounded solutions, and reconciles 
quite well the conflicting requirements of stability, 
accuracy and algorithmic simplicity. 

In what follows, we will present the details of 
the module and demonstrate how it is linked to the 
NPARC code. Application of the module will be 
reported in another paper 10 . 

2. Turbulence Models 

In accordance with the NPARC code, a non- 
dimensional form of equations is adopted. The three 
low Reynolds number two-equation turbulence mod- 
els built into the module have the following common 
form in which the Reynolds stresses t^(= —pu~Uj) 
are calculated by 


Tii = + Vi.i - I^m) - \pK6ii (2) 


2 



where R e is the reference Reynolds number. 

The 

2) Shih-Lumley’s K-e model. 


turbulent eddy viscosity /x t , the turbulent kinetic en- 
ergy K and its dissipation rate e are calculated by 
the following equations 

C„ = 0.09, Cj = 1.44, C 2 = 1.92 

(14) 

2 

Pt = R'fiiCpp— 

(3) 

o’K = 1, «r e = 1.3 



//i = [1 - exp (~aiR k - a s Rl - as-Rf)] 1 / 2 


{pK) t + [pUjK - R~ l {p + £-)Kj]j 

(4) 

R R'pK'fiy* 

P 

(15) 

= P- pe + D 


a x = 1.7 * 10" s , a s = 10" 9 , a 5 = 5 * 10~ 10 


(peh + \pu j <-R7 1 (p + %->Ai 


D = 0 

(16) 

U € 

€ € 2 
= fiCi—P — foCip— + E 

(5) 

E =??±SiSi 

pR 2 t ' ' 

(17) 

where 


S^sji 

(18) 

P = TijVij 

(6) 

= \(U XJ + U j>t ) - 

(19) 

f 1 = 1 

(7) 

The wall boundary conditions for K and e are 


h = 1 - 0.22exp[— (iZ t /6) 2 ] 

(8) 

K = 0.250u 2 

(20) 

lie 

(9) 

put 

€ = 0.251 

M 

(21) 


P is the turbulent production. For high Reynolds 
number models, / p = / x = f 2 = 1 and D = E = 0. 
The differences in the three models are given below: 


3) CMOTT K-e model. This model is the same as 
the Shih-Lumley model, except that is calculated 
by 


1) Chien’s X-e model. 

= 0.09, C x = 1.35, C 2 = 1.8 

( 10 ) 

<tk = 1, <r t - 1.3 

= 1 - exp(— 0.0115j/+) 

+ R e pu r y„ ( n ) 

y = 

p 

D = -2i?rV ~ (12) 

^ = _ 2 flrVe 

3/n 

At the wall, the values of iiT and e are both set to 
zero. 


= miii[0.09, {A 0 + A.U'K/e)- 1 ] (22) 

where 

Aq = 4, A, = \/6cos <{> (23) 


U' = (24) 

4> = iarccos(V6W0, W = (25) 

5* = ^, JVi = 5(Dij-^>) (26) 

In the above formulas, y* refers to the normal 
distance from the wall. 
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3. Calculation Procedure 


3.1 General Form of Turbulence Equations 

The turbulent transport equations ( 4) and ( 5) 
in general curvilinear coordinates (£, 77 ) may be writ- 
ten in the following form: 

d t {rJ- x p<t>) + d((rJ~ 1 pU<f>) + dr,{rJ- x pV<i > ) 

- 0 f (rJ- . V<fi) - d„(r/-V*V»7 • V<j>) (27) 
= rJ~ x S x 

where <f> stands for K or e, is the source terms on 
the right-hand side of Eqs.( 4) and ( 5), and 



(28) 

U = UU, + U 2 iy 

(29) 

V = U\T) X - 1 - 



(30) 

• V0 = Bl<f> ( + B\<t> v 

(31) 

Vr>-V<f>= B*<f> ( + B\ <f> v 


Bl = 


B\ = ZxVx + £yVy 

(32) 

B\ = T)x£ X + Tfy£y 


Bl = VxVz + rtyVy 



5| = d ( (rJ~ x B l 2 ^d v <f>) + d^rJ-'Blptdtf) (34) 

3.2 Numerical Discretization 

Equation ( 33) is discretized with the finite- 
volume approach. Integrating Eq.( 33) over a typical 
control volume centered at C (Fig.2) leads to a flux 
balance equation 

d t {p<f>)dil + F t -F w +F n -F, = J S+dQ (35) 

where Fi represents the total flux of <j> across the 
cell-face i(= e, w , 7i or s). Each of the surface fluxes 
Fi contains a convective contribution Ff and a dif- 
fusive contribution FP , that is 

Fi = Ff + F t P (36) 

Equation ( 35) involves no approximation and rep- 
resents a finite-volume analogue of the differential 
equation, Eq.( 33). 

Approximation of Convective Terms 

The convective contribution in Eq.( 36) can be 
approximated as 

Ff = Cifc (37) 

where Ci is the mass flux across the cell-face i, and 
can be calculated as 

C w = (rJ~ 1 pU) w 
C t = (rJ- x pU) e 

(38) 

C. = ( rJ~ x pV ), 


In the above equations, U\ and Ui are the velocity 
components along the x— and y— coordinates of the 
reference frame, r = 1 for planar and r = y for 
axisymmetric flows. 

Equation (27) can further be written as 

dt(rJ~ x p<f>) + d ( (r J~ x pU<P) + d v {rJ~ l pV<f>) 

- d((rJ~ x p.<fiBld(<l>) - d^TJ-’Lp.fBldr,#) (33) 

= rJ- x S x +Sl 

where Sj contains cross- derivative terms due to non- 
orthogonality 


C n = {rJ- l pV) n 

The calculation of fa is a key element for the 
accuracy and stability of numerical solutions. The 
more accurate schemes tend to be less stable, and 
vice versa. Here, the second-order accurate HLPA 
scheme 8 is used to calculate the face value of 4 > . This 
scheme is implemented in the following deferred cor- 
rection way proposed by Khosla and Rubin 11 

k = <f>? + \(4>f l - 1 -<l>*’ t - 1 ) (39) 

where u and h indicate the (first-order accurate) up- 
wind and higher-order schemes, 1—1 represents the 
previous iteration level, and A is a parameter which 
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blends the two schemes with limiting values A = 0 foi 
the upwind and A = 1 for the higher-order scheme. 
For the upwind scheme 


0£ = 4 U v <f>c 

<t>e = U+<j>c 4 U-<j> E 
= V+<i> s 4- V-<f>c 
« = W* 4 nrrfjr 


where 


D+rrl-DT, = 

^ + = i-vr, v;- = 

For the HLPA scheme 


0, if Ui > 0 

1, otherwise 


( 40 ) 


( 41 ) 


A <f>n = ^„ + a+(^ N - 4>c) 

+ K" ~ <t>N) 


<*+ = < 


1 if 


4c — 

— <£s 

4n — 4nN 

4>C ~ <f>NN 

4c - 2<j> w + 4>ww \ 


4>c - 4>ww 

0 otherwise 


1 if 


Ot.n ~ 


<f>W — 2<f>c + 4b 


4>w - <i>B 
0 otherwise 


< 1 


< 1 


I 1 * 

<f>B — 2 <j>c 4 4>w 

4>e — 4>w 

{ 0 otherwise 


[ 0, if Vi > 0 

(42) 


f 1 * 

i 

<f>c — 2 <f>E 4 <t>BE 

[ 1, otherwise 

= < 

<f>c - <t>EE 




0 otherwise 


< 1 


< 1 


( 47 ) 


( 48 ) 


( 49 ) 


4>% = 4w + U w 4>c + A<j> w 

4 h t = u+4> c + + A<fi' 


4) = v+<t>s + v~4>c + A <f>, 

<ft = V+<f>c + V n -<t> N + A<t> n 


( 43 ) 



r 

1 if 

i 

<t>c ~ %4>s 4 0ss 

a+ = < 

<t>C - 05S 


[ 0 otherwise 


< 1 



if 


4>s — 24c + 4>n 


4s — <t>N 

otherwise 


< 1 


( 50 ) 


where 


A<p w = U+a+(4>c — 4w) 
+ U~a~(4>w — 4c) 

A <f> t = U+a+(<j>E ~ <f>c) 
+ a~ (<f>c - 4 >e) 


4w — 4ww 
4c — 4>ww 
4c_ — 4b 
4>w - 4b 


$£ - 4w 
4 b — 4w 
4 b - 4eE 
4c — 4>bb 


( 44 ) 


( 45 ) 


A <f>, = V+a+(<f> c - 4s) 
+ V,~ a 7(4s - 4c) 


4s - 4s_S 
4c - 4ss 
4c - <I>N 
4s — 4 n 


( 46 ) 


1 if 


4 N — 2<t>C + 4s 


4n - 4s 

0 otherwise 


1 if 


4c — 2<j>s + 4nN 


a~ = < I 4c — 4 nn 

0 otherwise 


< 1 


< 1 


( 51 ) 


It can be seen that Eq.( 43) is in fact the result 
of the first-order upwinding with an additional term 
A 4>i added. The additional term may be viewed as 
an antidiffusive correction to the upwind scheme. 

The performance of the HLPA scheme is demon- 
strated in Figure 3 for a one-dimensional pure con- 
vection of a scalar profile. Figure 3(a) shows the test 
case: a rectangular S-profile at the initial time t=0 
is being transported in the x-direction by an invis- 
cid flow with prescribed velocities U=1 and V~0. 
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Figure 3(b) shows the computational domain. Cal- 
culations were carried out on two uniform grids, one 
comprising 201x2 and the other 1001x2 grid nodes. 
Two time increments were used, with A<=0.4 for the 
coarse and At=0.1 for the fine grid. Figures 3(c) 
and (d) show the S-profiles at the time t=100 calcu- 
lated by five different convection schemes: upwind, 
QUICK 12 , SOUCUP 13 , HLPA and SMART 14 . It is 
seen that HLPA captures steep gradients quite well 
while maintaining the boundedness of the solution, 
as compared to the other schemes. 

Approximation of Diffusive Terms 

The diffusive flux in Bq.( 36) can be divided into 
two parts 


FP = F t DN + FP° (52) 


The coefficients S£ and are defined such that 


S$> 0 

s$ < 0 


(56) 


This treatment enhances the stability of the numer- 
ical process and prevents the calculated value of <f> 
from becoming negative, especially in low turbulence 
regions. 

Approximation of Time Derivative 

Two time discretization schemes are used to ap- 
proximate the time derivative appearing in Eq.( 35). 
They are 

1) Fully Implicit (FI) scheme of first-order accu- 
racy which gives 


The first part FP n contains only one term which has 
the first derivative of <f> in the direction “normal” to 
the cell-face t. It can be written as 


F°" ■ =~D v (<t> c ~<t>w) 
f dn = _£> e ( 0 B - <f> c ) 

F ? N = -D,{4>c-4>s) 

F™ = -D n (4> N - <f>c) 

where 


d t {fxj>) = 


{p&F - {pf£_ 1 

At 


(57) 


2) Three-Level Fully Implicit (TLFI) scheme of 
second-order accuracy which gives 


Mp4>) = 


3 ( p 0 ) n - 4 ( p 0)"~ 1 + ( p ^)"- 2 
2 At 


( 58 ) 


where the superscript n refers to the current time 
level, and 


A t = t n - 1*- 1 


(59) 


Dtp — {t J 1 
D t = 

D $ = 

D n = (rJ- l Blti+) n 


( 54 ) 


The second part FP C involves the cross-derivative 
terms arising from grid non-orthogonality in 
Eq.( 34). Only the normal derivative diffusive flux, 
FP N y is coupled with the convective flux, Ff , to 
form the main coefficients of the difference equa- 
tion, while the cross-derivative diffusive flux, F ? c , 
is treated explicitly as a pseudo-source term to avoid 
the possibility of producing negative coefficients in 
an implicit treatment. 

Linearization of Source Terms 


Final Form of Discretized Equations 

After replacing all the terms in Eq.( 35) by their 
discretized analogues, we obtain the following differ- 
ence equation which relates the principal unknown 
<f>c to its four neighbours &(t = W y J E } 5, N) 


Ac<f>c = Ae<!>e + Aw<t>w + An4>s + As4>s 

+ Sp 


( 60 ) 


where 


Ab = D e — C e U e 
Aw = D„ + C„ V + 

A„ = D n - C n V~ (61) 


The source terms are linearized as 


As =D,+C,V+ 


S+ = S% +S$<f>c 


( 55 ) Ac = Ae + Aw + An + As — Sn 
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a. Initialize all field values. 


Sp = Sp + max(0, S£) + S| + S% 

(62) 

= Sjf + min(0, S|/<£) + Sj, 

where, Sp and Sjf aie the coefficients of the source 
term Sj, 

For the if-equation (<fi=K): 


Sj, = rJ _1 [max(0, P) + max(0, D)] 

( 63 ) 

S}f = rJ 1 [min(0, P) + min(0, D) — p4\/K 
For the e-equation (<f)=e): 


S}, = rJ -1 [max(0, fiCiP/K)e + max(0, E)\ 


S# = rJ -1 [min(0, hC^P/K) + min(0, E)/e ( 64 ) 

- fiCipe/K] 


S% is calculated by Eq.( 34), S| contains the deferred 
correction part in Eq.( 39) 


S% = -A [C e (# - tf) - - 4,1) 

+ C n (^-^“)-C' J (^-^)] 
Sj, and Sy are due to the use of 
FI scheme (Eq. 57): 


(65) 


c4 P jji— r 

- -^-tc 

C4 _ T*- l P 

bN ~ a r 


TLFI scheme (Eq. 58): 


(66) 


s% = r J^ 

p At 


( 2^' 1 - 0 . 5 ^“ 2 ) 


Sir = —1-5 


tJ l p 
At 


(67) 


In formulating Eq.( 60), the convection terms 
calculated by the upwind scheme are coupled with 
the normal diffusion terms to form the main coeffi- 
cients Ai. By this approach, the positivity of all the 
main coefficients is ensured so that the resulting co- 
efficient matrix will be always diagonally dominant. 


3.3 Calculation Sequence 

The sequence in which the calculation is carried 
out is as follows: 


b. Solve the U\- and ^-momentum equations using 

the NPARC code. 

c. Update the boundary conditions for the turbulent 

quantities. 

d. Solve the Jf-equation. 

e. Solve the e-equation. 

f. Update the turbulent eddy viscosity p t . 
g* Return to step b. 

The sequence of steps b to g are repeated until the 
convergence criterion is satisfied. The system of 
equations ( 60) are solved with the alternating di- 
rection TDMA of Thomas. 

4. Description of the Module 

To facilitate identification, all the subroutine 
names in the module start with CM. In order to 
use the module, the user only needs to call its two 
subroutines: CM AO and CMTURB, in the NPARC 
code. 

4.1 Subroutine CM AO 

The subroutine CMAO is for setting up the pa- 
rameters of the module. It contains the following 
arguments: 

JMAX,KMAX,NMAX,NTURB,IAXISY,RE, 

NJPAT,JPJ2,JPJM,JPK2,JPKM, 

NKPAT,KPJ2,KPJM,KPK2,KPKM, 

NJSEG,JLINE,JKLOW,JKHIGH,JTYPE, 

JSIGN,JEDGE, 

NKSEG, KLINE, KJLOW,KJHIGH,KTYPE, 
KSIGN,KEDGE 

These are all the variables already used in the 
NPARC code to define flow, geometric and bound- 
ary conditions. In addition, it has the following user- 
specified parameters: 

BDMAX(i), BDMIN(i) — Upper and lower bounds 
for the values of K (i=l), e (i=2) and (i=3). 
These bounds are introduced for numerical purposes 
only, that is, to prevent the corresponding turbu- 
lence quantities from becoming negative or abnor- 
mally large during the solution process. Currently, 
they are set to 

BDMAX(l)=1.0E-f6 

BDMAX(2)=1.0E+6 

BDMAX(3)=5.0E+3 

BDMIN(l)=1.0E-8 
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BDMIN(2)=1.0E-8 

BDMIN(3)=1.0E-3 

These should cover a wide range of the physically 
meaningful values of K, e and It is to be 

noted that these values are only valid for the non- 
dimension al turbulence quantities, as defined in the 
NPARC code. 

FDEFER — Blending factor A defined in Eq.( 39). 
Its value may vary from 0 to 1 with the limiting 
value 0 for the upwind and 1 for the ELPA scheme. 
The solution tends to be more stable, but also more 
diffusive when this factor is reduced. 

MODELS — Turbulence model selector with 1 for 
Chien, 2 for Shih-Lumley and 3 for CMOTT models. 
The Baldwin- Lomax model is used when NC is less 
than NTURB, which is defined in and is consistent 
with the NPARC code. 

RELAX(i) — Under-relaxation factors for K (i=l) 
and e (i=2). Currently, they are set to 
RELAX(1)=0.8 
RELAX(2)=0.8 

Should unstability occur, try to reduce these values. 

The subroutine CMAO needs to be called only 
once in the NPARC code. Currently, it is called in 
the subroutine INITIA as shown below 

SUBROUTINE INITIA 


DO 210 MB=l,NBLOCK 


CALL CMA0(..-) 
210 CONTINUE 


END 

4.2 Subroutine CMTURB 

This is the major interface between the NPARC 
and module. The user may simply replace all CALL 
MUTURB by CALL CMTURB in NAPRC. The 
subroutine MUTURB calls all turbulence model sub- 
routines in NPARC. 
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Fig. 2 A typical control volume centered at node C 
and its neighbouring nodes 
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